Mixed-Effects Model for Irradiance Time Above Ek
Angela Richards Donà
11/10/2022
#script for mixed model analysis of Ek and irradiance
Load the various libraries
library(lme4)## Loading required package: Matrixlibrary(lmerTest)##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## steplibrary(effects)## Loading required package: carData## lattice theme set by effectsTheme()
## See ?effectsTheme for details.library(car)
library(MuMIn)
library (dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:car':
##
## recode## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(emmeans)
library(DHARMa)## This is DHARMa 0.4.5. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')library(performance)
library(patchwork)
library(rstatix)##
## Attaching package: 'rstatix'## The following object is masked from 'package:stats':
##
## filter#for plots and tables
library(ggplot2)
library(ggpubr)
library(forcats)
library(RColorBrewer)
library(tidyverse)## ── Attaching packages
## ───────────────────────────────────────
## tidyverse 1.3.2 ──## ✔ tibble 3.1.8 ✔ purrr 0.3.5
## ✔ tidyr 1.2.1 ✔ stringr 1.4.1
## ✔ readr 2.1.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ rstatix::filter() masks dplyr::filter(), stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ tidyr::pack() masks Matrix::pack()
## ✖ dplyr::recode() masks car::recode()
## ✖ purrr::some() masks car::some()
## ✖ tidyr::unpack() masks Matrix::unpack()library(sjPlot)
library(sjmisc)##
## Attaching package: 'sjmisc'
##
## The following object is masked from 'package:purrr':
##
## is_empty
##
## The following object is masked from 'package:tidyr':
##
## replace_na
##
## The following object is masked from 'package:tibble':
##
## add_caselibrary(mmtable2)##
## Attaching package: 'mmtable2'
##
## The following object is masked from 'package:tidyr':
##
## table1library(gt)
library(purrr)
library(stringr)
library(tidyr)Load the data
ek_irrad_data <- read.csv("input_data/mean_supersaturation_by_period.csv")Assign run as a factor
ek_irrad_data$Run <- as.factor(ek_irrad_data$Run)Assign temperature as a factor
ek_irrad_data$Temperature <- as.factor(ek_irrad_data$Temp...C.)Assigns treatment as characters from integers then to factors
ek_irrad_data$Treatment <- as.factor(as.character(ek_irrad_data$Treatment))Toggle between the species for output. Use Day 9 for final analysis recent change: removing the very odd ek.1 value of 559.4 in hypnea dataset
hypnea <- subset(ek_irrad_data, Species == "hm")
hypnea$treatment_graph[hypnea$Treatment == 0] <- "1) 35ppt/0.5umol"
hypnea$treatment_graph[hypnea$Treatment == 1] <- "2) 35ppt/14umol"
hypnea$treatment_graph[hypnea$Treatment == 2] <- "3) 28ppt/27umol"
hypnea$treatment_graph[hypnea$Treatment == 3] <- "5) 18ppt/53umol"
hypnea$treatment_graph[hypnea$Treatment == 4] <- "6) 11ppt/80umol"
hypnea$treatment_graph[hypnea$Treatment == 2.5] <- "4) 28ppt/53umol"
ulva <- subset(ek_irrad_data, Species == "ul")
ulva$treatment_graph[ulva$Treatment == 0] <- "1) 35ppt/0.5umol"
ulva$treatment_graph[ulva$Treatment == 1] <- "2) 35ppt/14umol"
ulva$treatment_graph[ulva$Treatment == 2] <- "3) 28ppt/27umol"
ulva$treatment_graph[ulva$Treatment == 3] <- "5) 18ppt/53umol"
ulva$treatment_graph[ulva$Treatment == 4] <- "6) 11ppt/80umol"
ulva$treatment_graph[ulva$Treatment == 2.5] <- "4) 28ppt/53umol"Add growth rate from other dataset to this one and subset by species
growth_rate <- read.csv("/Users/Angela/src/work/limu/algal_growth_photosynthesis/data_input/all_runs_growth_102222.csv")
growth_rate$Species <- as.factor(growth_rate$Species)
growth_rate$treatment <- as.factor(growth_rate$treatment)Subset the growth rate data by Species
gr_ulva <- subset(growth_rate, Species == "Ul")
gr_hypnea <- subset(growth_rate, Species == "Hm" & final.weight != 0.1017)Calculate the growth rate from the data
ulva$growth_rate <- round((gr_ulva$final.weight - gr_ulva$Initial.weight) / gr_ulva$Initial.weight * 100, digits = 2)
hypnea$growth_rate <- round((gr_hypnea$final.weight - gr_hypnea$Initial.weight) / gr_hypnea$Initial.weight * 100, digits = 2)#_______________________ULVA____________________________
Run model without interaction between the treatments and temperature - irradiance_over_ek is in minutes. Take RLC.Order out of the random effects because it is causing problems of singularity. R2 is same with or without (+ (1 | RLC.Order))
ek_irrad_model_ulva <- lmer(formula = supersat_total ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID), data = ulva)Get an r2 for the model
r.squaredGLMM(ek_irrad_model_ulva)## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.## R2m R2c
## [1,] 0.3642221 0.8897065Make histograms of data and residuals (first test for mixed model use) and residual plots of the data for ulva
hist(ulva$supersat_total, main = paste("Ulva lactuca"), col = "olivedrab3", labels = TRUE)
hist(resid(ek_irrad_model_ulva))
plot(resid(ek_irrad_model_ulva) ~ fitted(ek_irrad_model_ulva))
qqnorm(resid(ek_irrad_model_ulva))
qqline(resid(ek_irrad_model_ulva))
ulva %>% ggplot(aes(supersat_total)) +
geom_histogram(binwidth=5, fill = "#5BB300", color = "black", size = 0.25, alpha = 0.85) +
theme_bw()
Check the performance of the model and run a summary
performance::check_model(ek_irrad_model_ulva)
summary(ek_irrad_model_ulva)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: supersat_total ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID)
## Data: ulva
##
## REML criterion at convergence: 2813
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1222 -0.4370 0.0455 0.4915 2.8439
##
## Random effects:
## Groups Name Variance Std.Dev.
## Plant.ID (Intercept) 655.3 25.60
## Run (Intercept) 2719.5 52.15
## Residual 708.3 26.61
## Number of obs: 288, groups: Plant.ID, 144; Run, 8
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 477.5252 37.4668 5.0114 12.745 5.21e-05 ***
## Treatment1 -50.3702 44.2995 4.9966 -1.137 0.3071
## Treatment2 -47.6295 44.2995 4.9966 -1.075 0.3315
## Treatment2.5 -157.9977 64.3119 4.8337 -2.457 0.0592 .
## Treatment3 -80.6781 44.2995 4.9966 -1.821 0.1283
## Treatment4 -88.4096 44.2995 4.9966 -1.996 0.1025
## Temperature27 0.1846 6.8448 115.6816 0.027 0.9785
## Temperature30 -2.3436 6.8448 115.6816 -0.342 0.7327
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trtmn1 Trtmn2 Trt2.5 Trtmn3 Trtmn4 Tmpr27
## Treatment1 -0.836
## Treatment2 -0.836 0.992
## Treatmnt2.5 -0.576 0.487 0.487
## Treatment3 -0.836 0.992 0.992 0.487
## Treatment4 -0.836 0.992 0.992 0.487 0.992
## Temperatr27 -0.091 0.000 0.000 0.000 0.000 0.000
## Temperatr30 -0.091 0.000 0.000 0.000 0.000 0.000 0.500Check for equal variance to determine which ANOVA to use
bartlett.test(supersat_total ~ Treatment, data = ulva)##
## Bartlett test of homogeneity of variances
##
## data: supersat_total by Treatment
## Bartlett's K-squared = 19.346, df = 5, p-value = 0.001656Run Welch’s ANOVA if not equal variance. Code forces one independent variable at a time. Games Howell Test for pairwise comparison for the only independent variable that was found to have significant results.
welch_anova_treatment_ulva <- oneway.test(supersat_total ~ Treatment, data = ulva, var.equal = FALSE)
welch_anova_treatment_ulva##
## One-way analysis of means (not assuming equal variances)
##
## data: supersat_total and Treatment
## F = 74.182, num df = 5.00, denom df = 130.53, p-value < 2.2e-16welch_anova_temp_ulva <- oneway.test(supersat_total ~ Temperature, data = ulva, var.equal = FALSE)
welch_anova_temp_ulva##
## One-way analysis of means (not assuming equal variances)
##
## data: supersat_total and Temperature
## F = 0.16463, num df = 2.00, denom df = 189.68, p-value = 0.8483games_howell_test(ulva, supersat_total ~ Treatment, conf.level = 0.95, detailed = FALSE)## # A tibble: 15 × 8
## .y. group1 group2 estimate conf.low conf.high p.adj p.adj.sig…¹
## * <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 supersat_total 0 1 -44.9 -77.4 -12.4 2 e- 3 **
## 2 supersat_total 0 2 -42.1 -74.8 -9.53 4 e- 3 **
## 3 supersat_total 0 2.5 -158. -183. -133. 4.75e-10 ****
## 4 supersat_total 0 3 -75.2 -108. -42.2 4.24e- 8 ****
## 5 supersat_total 0 4 -82.9 -115. -50.4 1.34e- 9 ****
## 6 supersat_total 1 2 2.74 -34.4 39.9 1 e+ 0 ns
## 7 supersat_total 1 2.5 -113. -144. -82.2 0 ****
## 8 supersat_total 1 3 -30.3 -67.8 7.23 1.85e- 1 ns
## 9 supersat_total 1 4 -38.0 -75.1 -0.939 4.1 e- 2 *
## 10 supersat_total 2 2.5 -116. -147. -84.8 0 ****
## 11 supersat_total 2 3 -33.0 -70.7 4.59 1.19e- 1 ns
## 12 supersat_total 2 4 -40.8 -78.0 -3.58 2.3 e- 2 *
## 13 supersat_total 2.5 3 82.8 51.3 114. 5.32e-10 ****
## 14 supersat_total 2.5 4 75.1 44.1 106. 7.93e- 9 ****
## 15 supersat_total 3 4 -7.73 -45.3 29.8 9.91e- 1 ns
## # … with abbreviated variable name ¹p.adj.signifplot(allEffects(ek_irrad_model_ulva))
Summarize the means for rETRmax
ulva %>% group_by(Treatment) %>% summarise_at(vars(supersat_total), list(mean = mean))## # A tibble: 6 × 2
## Treatment mean
## <fct> <dbl>
## 1 0 477.
## 2 1 432.
## 3 2 435.
## 4 2.5 319.
## 5 3 402.
## 6 4 394.ulva %>% group_by(Run) %>% summarise_at(vars(supersat_total), list(mean = mean))## # A tibble: 8 × 2
## Run mean
## <fct> <dbl>
## 1 1 461.
## 2 2 465.
## 3 3 422.
## 4 3.5 352.
## 5 4 349.
## 6 6 452.
## 7 6.5 501.
## 8 9 319.Plot a regression between the photosynthetic independent variables of interest and growth rate
ulva_growth_irrad_ek_graph <- ggplot(ulva, aes(x=supersat_total, y=growth_rate)) +
geom_point(alpha = 0.5, size = 3, aes(color = Treatment), show.legend = TRUE) +
geom_smooth(method = "lm", col = "black") + theme_bw() +
labs(title = "A", subtitle = "Ulva lactuca", x = "Mean saturation time (min)",
y = "growth rate (%)") + stat_regline_equation(label.x = 450, label.y = 225) +
stat_cor(label.x = 450, label.y = 215) +
theme(plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05),
plot.subtitle = element_text(face = "italic", vjust = -20, hjust = 0.05))
ulva_growth_irrad_ek_graph## `geom_smooth()` using formula 'y ~ x'
Plot a histogram in ggplot2
ulva %>% ggplot(aes(supersat_total)) +
geom_histogram(binwidth=5, fill = "#5BB300", color = "black", size = 0.25, alpha = 0.85) +
theme_bw()
Make a boxplot for publications of the results
ulva %>% ggplot(aes(treatment_graph, supersat_total)) +
geom_boxplot(size=0.5) +
geom_point(alpha = 0.5, size = 3, aes(color = Temperature), show.legend = FALSE) +
labs(x="Treatment (salinty/nitrate)", y= "Mean Daily Saturation Time (min)", title= "A", subtitle = "Ulva lactuca") +
scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "28ppt/53umolN", "18ppt/53umolN", "11ppt/80umolN")) +
ylim(100, 600) + stat_mean() +
geom_hline(yintercept=100, color = "orange", size = 0.5, alpha = 0.5) +
scale_color_manual(values = c("#999999", "#E69F00", "#56B4E9")) +
theme_bw() +
theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05),
plot.subtitle = element_text(face = "italic", vjust = -20, hjust = 0.05))
##Run model for relative Hsat
rel_hsat_model_ulva <- lmer(formula = supersat_rel ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID), data = ulva)hist(ulva$supersat_rel, main = paste("Ulva lactuca"), col = "olivedrab3", labels = TRUE)
hist(resid(rel_hsat_model_ulva))
plot(resid(rel_hsat_model_ulva) ~ fitted(rel_hsat_model_ulva))
qqnorm(resid(rel_hsat_model_ulva))
qqline(resid(rel_hsat_model_ulva))
ulva %>% ggplot(aes(supersat_total)) +
geom_histogram(binwidth=5, fill = "#5BB300", color = "black", size = 0.25, alpha = 0.85) +
theme_bw()
Check the performance of the model
performance::check_model(rel_hsat_model_ulva)
r.squaredGLMM(rel_hsat_model_ulva)## R2m R2c
## [1,] 0.4420469 0.8743361summary(rel_hsat_model_ulva)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: supersat_rel ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID)
## Data: ulva
##
## REML criterion at convergence: -802
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2462 -0.4333 0.0519 0.4954 2.7114
##
## Random effects:
## Groups Name Variance Std.Dev.
## Plant.ID (Intercept) 0.001543 0.03928
## Run (Intercept) 0.004649 0.06818
## Residual 0.001800 0.04242
## Number of obs: 288, groups: Plant.ID, 144; Run, 8
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.751352 0.049313 5.008835 15.236 2.18e-05 ***
## Treatment1 -0.083563 0.058292 4.988804 -1.434 0.2113
## Treatment2 -0.079353 0.058292 4.988804 -1.361 0.2317
## Treatment2.5 -0.260459 0.084333 4.760091 -3.088 0.0290 *
## Treatment3 -0.131558 0.058292 4.988804 -2.257 0.0738 .
## Treatment4 -0.143701 0.058292 4.988804 -2.465 0.0570 .
## Temperature27 0.009256 0.010659 112.449439 0.868 0.3870
## Temperature30 0.010096 0.010659 112.449439 0.947 0.3456
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trtmn1 Trtmn2 Trt2.5 Trtmn3 Trtmn4 Tmpr27
## Treatment1 -0.833
## Treatment2 -0.833 0.989
## Treatmnt2.5 -0.576 0.487 0.487
## Treatment3 -0.833 0.989 0.989 0.487
## Treatment4 -0.833 0.989 0.989 0.487 0.989
## Temperatr27 -0.108 0.000 0.000 0.000 0.000 0.000
## Temperatr30 -0.108 0.000 0.000 0.000 0.000 0.000 0.500Check for equal variance
bartlett.test(supersat_rel ~ Treatment, data = ulva)##
## Bartlett test of homogeneity of variances
##
## data: supersat_rel by Treatment
## Bartlett's K-squared = 10.634, df = 5, p-value = 0.05915#run Welch's ANOVA if not equal variance
welch_anova_treatment_ulva <- oneway.test(supersat_rel ~ Treatment, data = ulva, var.equal = FALSE)
welch_anova_treatment_ulva##
## One-way analysis of means (not assuming equal variances)
##
## data: supersat_rel and Treatment
## F = 84.949, num df = 5.00, denom df = 131.02, p-value < 2.2e-16welch_anova_temp_ulva <- oneway.test(supersat_rel ~ Temperature, data = ulva, var.equal = FALSE)
welch_anova_temp_ulva##
## One-way analysis of means (not assuming equal variances)
##
## data: supersat_rel and Temperature
## F = 0.5808, num df = 2.00, denom df = 189.75, p-value = 0.5604games_howell_test(ulva, supersat_rel ~ Treatment, conf.level = 0.95, detailed = TRUE)## # A tibble: 15 × 14
## .y. group1 group2 n1 n2 estimate conf.…¹ conf.h…² se stati…³
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 supersat… 0 1 48 48 -0.0770 -0.124 -3.06e-2 0.0113 4.83
## 2 supersat… 0 2 48 48 -0.0728 -0.119 -2.63e-2 0.0113 4.56
## 3 supersat… 0 2.5 48 48 -0.260 -0.299 -2.22e-1 0.00943 19.5
## 4 supersat… 0 3 48 48 -0.125 -0.173 -7.72e-2 0.0116 7.62
## 5 supersat… 0 4 48 48 -0.137 -0.184 -9.01e-2 0.0114 8.49
## 6 supersat… 1 2 48 48 0.00421 -0.0463 5.47e-2 0.0123 0.243
## 7 supersat… 1 2.5 48 48 -0.183 -0.227 -1.40e-1 0.0106 12.3
## 8 supersat… 1 3 48 48 -0.0480 -0.0996 3.66e-3 0.0126 2.70
## 9 supersat… 1 4 48 48 -0.0601 -0.111 -9.15e-3 0.0124 3.43
## 10 supersat… 2 2.5 48 48 -0.188 -0.231 -1.44e-1 0.0106 12.5
## 11 supersat… 2 3 48 48 -0.0522 -0.104 -4.81e-4 0.0126 2.94
## 12 supersat… 2 4 48 48 -0.0643 -0.115 -1.33e-2 0.0124 3.67
## 13 supersat… 2.5 3 48 48 0.135 0.0904 1.80e-1 0.0109 8.77
## 14 supersat… 2.5 4 48 48 0.123 0.0790 1.68e-1 0.0107 8.12
## 15 supersat… 3 4 48 48 -0.0121 -0.0643 4.00e-2 0.0127 0.677
## # … with 4 more variables: df <dbl>, p.adj <dbl>, p.adj.signif <chr>,
## # method <chr>, and abbreviated variable names ¹conf.low, ²conf.high,
## # ³statisticplot(allEffects(rel_hsat_model_ulva))
ulva %>% ggplot(aes(treatment_graph, supersat_rel)) +
geom_boxplot(size=0.5) +
geom_point(alpha = 0.5, size = 3, aes(color = Temperature), show.legend = FALSE) +
labs(x="Treatment (salinty/nitrate)", y= "Mean Rel. Hsat Time", title= "A", subtitle = "Ulva lactuca") +
scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "28ppt/53umolN", "18ppt/53umolN", "11ppt/80umolN")) +
ylim(0, 1) + stat_mean() +
geom_hline(yintercept=0.5, color = "orange", size = 0.5, alpha = 0.5) +
scale_color_manual(values = c("#999999", "#E69F00", "#56B4E9")) +
theme_bw() +
theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05),
plot.subtitle = element_text(face = "italic", vjust = -20, hjust = 0.05))
Summarize the means for relative Hsat
ulva %>% group_by(Treatment) %>% summarise_at(vars(supersat_rel), list(mean = mean))## # A tibble: 6 × 2
## Treatment mean
## <fct> <dbl>
## 1 0 0.758
## 2 1 0.681
## 3 2 0.685
## 4 2.5 0.497
## 5 3 0.633
## 6 4 0.621ulva %>% group_by(Run) %>% summarise_at(vars(supersat_rel), list(mean = mean))## # A tibble: 8 × 2
## Run mean
## <fct> <dbl>
## 1 1 0.704
## 2 2 0.723
## 3 3 0.674
## 4 3.5 0.567
## 5 4 0.572
## 6 6 0.723
## 7 6.5 0.792
## 8 9 0.497ulva %>% group_by(Run) %>% summarise_at(vars(day_length_avg), list(mean = mean))## # A tibble: 8 × 2
## Run mean
## <fct> <dbl>
## 1 1 654.
## 2 2 643.
## 3 3 626.
## 4 3.5 621.
## 5 4 611.
## 6 6 626.
## 7 6.5 633.
## 8 9 641.ulva_growth_rel_hsat_graph <- ggplot(ulva, aes(x=supersat_rel, y=growth_rate)) +
geom_point(alpha = 0.5, size = 3, aes(color = Treatment), show.legend = TRUE) +
geom_smooth(method = "lm", col = "black") + theme_bw() +
labs(title = "A", subtitle = "Ulva lactuca", x = "Mean Rel. Hsat",
y = "growth rate (%)") + stat_regline_equation(label.x = 0.68, label.y = 200) +
stat_cor(label.x = 0.68, label.y = 215) +
theme(plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05),
plot.subtitle = element_text(face = "italic", vjust = -20, hjust = 0.05))
ulva_growth_rel_hsat_graph## `geom_smooth()` using formula 'y ~ x'
#______________________HYPNEA___________________________
Run model for Hypnea
ek_irrad_model_hypnea <- lmer(formula = supersat_total ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID), data = hypnea)Get an r2 for the model
r.squaredGLMM(ek_irrad_model_hypnea)## R2m R2c
## [1,] 0.1460739 0.7597093Make histograms of data and residuals (first test for mixed model use) and residual plots of the data for Hypnea
hist(hypnea$supersat_total, main = paste("Hypnea musciformis"), col = "maroon", labels = TRUE)
hist(resid(ek_irrad_model_hypnea))
plot(resid(ek_irrad_model_hypnea) ~ fitted(ek_irrad_model_hypnea))
qqnorm(resid(ek_irrad_model_hypnea))
qqline(resid(ek_irrad_model_hypnea))
Plot a histogram in ggplot
hypnea %>% ggplot(aes(supersat_total)) +
geom_histogram(binwidth=5, fill = "#a8325e", color = "black", size = 0.25, alpha = 0.85) +
theme_bw()
Check the performance of the model
performance::check_model(ek_irrad_model_hypnea)
Run a summary of the model
summary(ek_irrad_model_hypnea)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: supersat_total ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID)
## Data: hypnea
##
## REML criterion at convergence: 2904.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3115 -0.4540 0.0553 0.5465 2.1801
##
## Random effects:
## Groups Name Variance Std.Dev.
## Plant.ID (Intercept) 787.4 28.06
## Run (Intercept) 2032.5 45.08
## Residual 1104.2 33.23
## Number of obs: 287, groups: Plant.ID, 144; Run, 8
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 396.535 32.811 5.043 12.086 6.47e-05 ***
## Treatment1 -49.539 38.782 5.021 -1.277 0.257345
## Treatment2 -35.766 38.782 5.021 -0.922 0.398558
## Treatment2.5 -37.009 55.925 4.729 -0.662 0.538991
## Treatment3 -51.386 38.782 5.021 -1.325 0.242263
## Treatment4 -69.747 38.789 5.025 -1.798 0.131789
## Temperature27 29.394 7.910 131.267 3.716 0.000299 ***
## Temperature30 32.267 7.915 131.714 4.077 7.86e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trtmn1 Trtmn2 Trt2.5 Trtmn3 Trtmn4 Tmpr27
## Treatment1 -0.830
## Treatment2 -0.830 0.985
## Treatmnt2.5 -0.575 0.487 0.487
## Treatment3 -0.830 0.985 0.985 0.487
## Treatment4 -0.830 0.985 0.985 0.487 0.985
## Temperatr27 -0.121 0.000 0.000 0.000 0.000 0.000
## Temperatr30 -0.121 0.000 0.000 0.000 0.000 0.001 0.500Check for equal variance
bartlett.test(supersat_total ~ Treatment, data = hypnea)##
## Bartlett test of homogeneity of variances
##
## data: supersat_total by Treatment
## Bartlett's K-squared = 50.954, df = 5, p-value = 8.84e-10Run Welch’s ANOVA if not equal variance. Code forces one independent variable at a time. Games Howell Test for pairwise comparison for the only independent variable that was found to have significant results.
welch_anova_treatment <- oneway.test(supersat_total ~ Treatment, data = hypnea, var.equal = FALSE)
welch_anova_treatment##
## One-way analysis of means (not assuming equal variances)
##
## data: supersat_total and Treatment
## F = 7.4774, num df = 5.00, denom df = 127.19, p-value = 3.509e-06welch_anova_temp <- oneway.test(supersat_total ~ Temperature, data = hypnea, var.equal = FALSE)
welch_anova_temp##
## One-way analysis of means (not assuming equal variances)
##
## data: supersat_total and Temperature
## F = 11.458, num df = 2.00, denom df = 184.95, p-value = 2.036e-05games_howell_test(hypnea, supersat_total ~ Treatment, conf.level = 0.95, detailed = FALSE)## # A tibble: 15 × 8
## .y. group1 group2 estimate conf.low conf.high p.adj p.adj.s…¹
## * <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 supersat_total 0 1 -42.7 -81.9 -3.53 0.025 *
## 2 supersat_total 0 2 -28.9 -63.5 5.60 0.153 ns
## 3 supersat_total 0 2.5 -37.0 -60.0 -14.0 0.000153 ***
## 4 supersat_total 0 3 -44.6 -79.3 -9.82 0.004 **
## 5 supersat_total 0 4 -62.0 -94.5 -29.6 0.00000416 ****
## 6 supersat_total 1 2 13.8 -30.4 57.9 0.944 ns
## 7 supersat_total 1 2.5 5.71 -30.6 42.0 0.997 ns
## 8 supersat_total 1 3 -1.85 -46.2 42.5 1 ns
## 9 supersat_total 1 4 -19.3 -61.9 23.3 0.773 ns
## 10 supersat_total 2 2.5 -8.07 -39.2 23.1 0.973 ns
## 11 supersat_total 2 3 -15.6 -56.0 24.7 0.869 ns
## 12 supersat_total 2 4 -33.1 -71.5 5.37 0.134 ns
## 13 supersat_total 2.5 3 -7.55 -39.0 23.8 0.98 ns
## 14 supersat_total 2.5 4 -25.0 -53.8 3.79 0.125 ns
## 15 supersat_total 3 4 -17.5 -56.1 21.2 0.776 ns
## # … with abbreviated variable name ¹p.adj.signifgames_howell_test(hypnea, supersat_total ~ Temperature, conf.level = 0.95, detailed = FALSE)## # A tibble: 3 × 8
## .y. group1 group2 estimate conf.low conf.high p.adj p.adj.sig…¹
## * <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 supersat_total 20 27 38.3 16.4 60.2 0.000161 ***
## 2 supersat_total 20 30 43.3 21.0 65.6 0.0000254 ****
## 3 supersat_total 27 30 4.96 -12.7 22.6 0.784 ns
## # … with abbreviated variable name ¹p.adj.signifIf ANOVA more appropriate (not this case) run from code above with Tukey’s HSD pairwise comparisons where appropriate
Plot results for basic output
plot(allEffects(ek_irrad_model_hypnea))
Plot a regression between the photosynthetic independent variables of interest and growth rate and a histogram of the data
hypnea_growth_irrad_ek_graph <- ggplot(hypnea, aes(x=supersat_total, y=growth_rate)) +
geom_point(alpha = 0.5, size = 3, aes(color = Treatment), show.legend = TRUE) +
geom_smooth(method = "lm", col = "black") + theme_bw() +
labs(title = "B", subtitle = "Hypnea musciformis", x = "Mean saturation time (min)", y = "growth rate (%)") +
stat_regline_equation(label.x = 400, label.y = 150) + stat_cor(label.x = 400, label.y = 140) +
theme(plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05),
plot.subtitle = element_text(face = "italic", vjust = -20, hjust = 0.05))
hypnea_growth_irrad_ek_graph## `geom_smooth()` using formula 'y ~ x'
Make a boxplot for publications of the results
hypnea %>% ggplot(aes(treatment_graph, supersat_total)) +
geom_boxplot(size=0.5) +
geom_point(alpha = 0.5, size = 3, aes(color = Temperature), show.legend = TRUE) +
labs(x="Treatment (salinty/nitrate)", y= "Mean Daily Saturation Time (min)", title= "B", subtitle = "Hypnea musciformis") +
scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "28ppt/53umolN", "18ppt/53umolN", "11ppt/80umolN")) +
ylim(100, 600) + stat_mean() +
geom_hline(yintercept=100, color = "orange", size = 0.5, alpha = 0.5) +
scale_color_manual(values = c("#999999", "#E69F00", "#56B4E9")) +
theme_bw() +
theme(legend.position = c(0.88,0.88), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05),
plot.subtitle = element_text(face = "italic", vjust = -20, hjust = 0.05))
Summarize the means for rETRmax
hypnea %>% group_by(Treatment) %>% summarise_at(vars(supersat_total), list(mean = mean))## # A tibble: 6 × 2
## Treatment mean
## <fct> <dbl>
## 1 0 417.
## 2 1 374.
## 3 2 388.
## 4 2.5 380.
## 5 3 373.
## 6 4 355.hypnea %>% group_by(Run) %>% summarise_at(vars(supersat_total), list(mean = mean))## # A tibble: 8 × 2
## Run mean
## <fct> <dbl>
## 1 1 407.
## 2 2 426.
## 3 3 360.
## 4 3.5 312.
## 5 4 320.
## 6 7 380.
## 7 8 406.
## 8 9 428.##Run the model for Hsat
rel_hsat_model_hypnea <- lmer(formula = supersat_rel ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID), data = hypnea)Make a histogram and residual plots of the data for ulva
hist(hypnea$supersat_rel, main = paste("Hypnea musciformis"), col = "maroon", labels = TRUE)
hist(resid(rel_hsat_model_hypnea))
plot(resid(rel_hsat_model_hypnea) ~ fitted(rel_hsat_model_hypnea))
qqnorm(resid(rel_hsat_model_hypnea))
qqline(resid(rel_hsat_model_hypnea))
hypnea %>% ggplot(aes(supersat_total)) +
geom_histogram(binwidth=5, fill = "#a8325e", color = "black", size = 0.25, alpha = 0.85) +
theme_bw()
Check the performance of the model
performance::check_model(rel_hsat_model_hypnea)
r.squaredGLMM(rel_hsat_model_hypnea)## R2m R2c
## [1,] 0.2080092 0.7254956summary(rel_hsat_model_hypnea)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: supersat_rel ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID)
## Data: hypnea
##
## REML criterion at convergence: -701.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4632 -0.4603 0.0425 0.5323 2.3552
##
## Random effects:
## Groups Name Variance Std.Dev.
## Plant.ID (Intercept) 0.001822 0.04269
## Run (Intercept) 0.003371 0.05806
## Residual 0.002755 0.05249
## Number of obs: 287, groups: Plant.ID, 144; Run, 8
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.61738 0.04279 5.05849 14.429 2.65e-05 ***
## Treatment1 -0.07498 0.05056 5.02872 -1.483 0.1978
## Treatment2 -0.05364 0.05056 5.02872 -1.061 0.3369
## Treatment2.5 -0.10415 0.07243 4.61640 -1.438 0.2147
## Treatment3 -0.07807 0.05056 5.02872 -1.544 0.1828
## Treatment4 -0.10692 0.05057 5.03413 -2.114 0.0878 .
## Temperature27 0.05385 0.01223 133.00960 4.402 2.18e-05 ***
## Temperature30 0.06321 0.01224 133.48716 5.164 8.60e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trtmn1 Trtmn2 Trt2.5 Trtmn3 Trtmn4 Tmpr27
## Treatment1 -0.823
## Treatment2 -0.823 0.978
## Treatmnt2.5 -0.575 0.486 0.486
## Treatment3 -0.823 0.978 0.978 0.486
## Treatment4 -0.823 0.977 0.977 0.486 0.977
## Temperatr27 -0.143 0.000 0.000 0.000 0.000 0.000
## Temperatr30 -0.143 0.000 0.000 0.000 0.000 0.001 0.500Check for equal variance
bartlett.test(supersat_rel ~ Treatment, data = hypnea)##
## Bartlett test of homogeneity of variances
##
## data: supersat_rel by Treatment
## Bartlett's K-squared = 49.512, df = 5, p-value = 1.744e-09#run Welch's ANOVA if not equal variance
welch_anova_treatment <- oneway.test(supersat_rel ~ Treatment, data = hypnea, var.equal = FALSE)
welch_anova_treatment##
## One-way analysis of means (not assuming equal variances)
##
## data: supersat_rel and Treatment
## F = 16.025, num df = 5.00, denom df = 127.18, p-value = 3.143e-12welch_anova_temp <- oneway.test(supersat_rel ~ Temperature, data = hypnea, var.equal = FALSE)
welch_anova_temp##
## One-way analysis of means (not assuming equal variances)
##
## data: supersat_rel and Temperature
## F = 16.829, num df = 2.00, denom df = 185.09, p-value = 1.927e-07games_howell_test(hypnea, supersat_rel ~ Treatment, conf.level = 0.95, detailed = TRUE)## # A tibble: 15 × 14
## .y. group1 group2 n1 n2 estimate conf.…¹ conf.h…² se stati…³
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 supersat… 0 1 48 48 -0.0663 -0.126 -0.00664 0.0144 3.25
## 2 supersat… 0 2 48 48 -0.0449 -0.0959 0.00601 0.0124 2.57
## 3 supersat… 0 2.5 48 48 -0.104 -0.140 -0.0687 0.00858 8.58
## 4 supersat… 0 3 48 48 -0.0694 -0.122 -0.0171 0.0127 3.87
## 5 supersat… 0 4 48 47 -0.0973 -0.146 -0.0488 0.0118 5.85
## 6 supersat… 1 2 48 48 0.0213 -0.0440 0.0866 0.0159 0.952
## 7 supersat… 1 2.5 48 48 -0.0379 -0.0926 0.0168 0.0131 2.04
## 8 supersat… 1 3 48 48 -0.00309 -0.0694 0.0632 0.0161 0.136
## 9 supersat… 1 4 48 47 -0.0310 -0.0944 0.0324 0.0154 1.42
## 10 supersat… 2 2.5 48 48 -0.0592 -0.104 -0.0143 0.0108 3.87
## 11 supersat… 2 3 48 48 -0.0244 -0.0832 0.0343 0.0143 1.21
## 12 supersat… 2 4 48 47 -0.0523 -0.108 0.00310 0.0135 2.75
## 13 supersat… 2.5 3 48 48 0.0348 -0.0117 0.0812 0.0112 2.20
## 14 supersat… 2.5 4 48 47 0.00689 -0.0351 0.0489 0.0101 0.482
## 15 supersat… 3 4 48 47 -0.0279 -0.0845 0.0287 0.0138 1.43
## # … with 4 more variables: df <dbl>, p.adj <dbl>, p.adj.signif <chr>,
## # method <chr>, and abbreviated variable names ¹conf.low, ²conf.high,
## # ³statisticgames_howell_test(hypnea, supersat_rel ~ Temperature, conf.level = 0.95, detailed = TRUE)## # A tibble: 3 × 14
## .y. group1 group2 n1 n2 estim…¹ conf.…² conf.…³ se stati…⁴ df
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 super… 20 27 96 96 0.0674 0.0350 0.0999 0.00971 4.91 166.
## 2 super… 20 30 96 95 0.0785 0.0455 0.112 0.00989 5.62 171.
## 3 super… 27 30 96 95 0.0111 -0.0152 0.0374 0.00786 0.998 188.
## # … with 3 more variables: p.adj <dbl>, p.adj.signif <chr>, method <chr>, and
## # abbreviated variable names ¹estimate, ²conf.low, ³conf.high, ⁴statisticplot(allEffects(rel_hsat_model_hypnea))
hypnea %>% ggplot(aes(treatment_graph, supersat_rel)) +
geom_boxplot(size=0.5) +
geom_point(alpha = 0.5, size = 3, aes(color = Temperature), show.legend = TRUE) +
labs(x="Treatment (salinty/nitrate)", y= "Mean Rel. Hsat Time", title= "B", subtitle = "Hypnea musciformis") +
scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "28ppt/53umolN", "18ppt/53umolN", "11ppt/80umolN")) +
ylim(0, 1) + stat_mean() +
geom_hline(yintercept=0.5, color = "orange", size = 0.5, alpha = 0.5) +
scale_color_manual(values = c("#999999", "#E69F00", "#56B4E9")) +
theme_bw() +
theme(legend.position = c(0.88,0.88), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05),
plot.subtitle = element_text(face = "italic", vjust = -20, hjust = 0.05))
hypnea_growth_rel_hsat_graph <- ggplot(hypnea, aes(x=supersat_rel, y=growth_rate)) +
geom_point(alpha = 0.5, size = 3, aes(color = Treatment), show.legend = TRUE) +
geom_smooth(method = "lm", col = "black") + theme_bw() +
labs(title = "B", subtitle = "Hypnea musciformis", x = "Mean Hsat Time (min)", y = "growth rate (%)") +
stat_regline_equation(label.x = 0.68, label.y = 100) + stat_cor(label.x = 0.68, label.y = 108) +
theme(plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05),
plot.subtitle = element_text(face = "italic", vjust = -20, hjust = 0.05))
hypnea_growth_rel_hsat_graph## `geom_smooth()` using formula 'y ~ x'